#include <stdio.h>
#include <stdlib.h>
#include <math.h>
// #include <atlas/clapack.h>
#include "clapack.h"
#include "spline.h"


/*define and create the cubic spline between lowergap and uppergap*/
void splinefunc(double spline[] ,int length, double zi, double zj, double delta_z){

    int i;
    //printf(" %f", pow(zi,0.0));
   /* double mat_reg[4*4] =  { 1.0 ,zi, pow(zi,2.0), pow(zi,3.0),
                             1.0, zj, pow(zj,2.0), pow(zj,3.0),
                             0.0, 1.0, 2.0*zi, 3.0*pow(zi,2.0),
                             0.0, 1.0, 2.0*zj, 3.0*pow(zj,2.0)};*/
     double mat_reg[4*4] =  { 1.0 , 1.0,0.0,0.0,
                              zi,zj,1.0,1.0,
                              pow(zi,2.0),pow(zj,2.0),2.0*zi,2.0*zj,
                              pow(zi,3.0) , pow(zj,3.0),3*pow(zi,2.0),3*pow(zj,2.0)};

    double b[4]          =  {1.0 , 0.0, 0.0 , 0.0};



    //res = mat_reg\b;
    integer n = 4;
    integer nhrs = 1;
    integer lda = 4;
    integer ldb = 4;
    integer ipiv[4];
    integer info;
    dgesv_(&n, &nhrs, mat_reg, &lda, ipiv, b, &ldb, &info);


    for(i = 0; i<length; i++){
        //spline[i] = b[1]+b[2]*(zi+i*delta_z)+b[3]*pow(zi+i*delta_z,2.0)+b[4]*pow(zi+i*delta_z,3.0);
        spline[i] = b[0]+b[1]*(zi+(i+1)*delta_z)+b[2]*pow(zi+(i+1)*delta_z,2.0)+b[3]*pow(zi+(i+1)*delta_z,3.0);
    }
}
